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We propose a method to extract from experimental data the subdiffusion parameter a and subdif- 
fusion coefficient D a which are defined by means of the relation (x 2 ) — r ^_° ) t a where {x 2 } denotes 
a mean square displacement of a random walker starting from x = at the initial time t = 0. The 
method exploits a membrane system where a substance of interest is transported in a solvent from 
QQ ' one vessel to another across a thin membrane which plays here only an auxiliary role. Using such a 

£SJ ' system, we experimentally study a diffusion of glucose and sucrose in a gel solvent. We find a fully 

analytic solution of the fractional subdiffusion equation with the initial and boundary conditions 
representing the system under study. Confronting the experimental data with the derived formulas, 
we show a subdiffusive character of the sugar transport in gel solvent. We precisely determine the 
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. parameter a, which is smaller than 1, and the subdiffusion coefficient D a . 

I. INTRODUCTION 

•*-> 

Subdiffusion occurs in various systems. As examples, we mention here a diffusion in porous media or charge carriers 
transport in amorphous semiconductors [1-3]. The subdiffusion is characterized by a time dependence of the mean 
square displacement of a Brownian particle. When the particle starts form x = at the initial time t = this 
I , dependence in a one-dimension system is of the form 
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where D a is the subdiffusion coefficient measured in the units [m 2 /s a ]; the parameter a, which we call here as a 
i subdiffusion parameter, obeys < a < 1. For a = 1 one deals with the normal or Gaussian diffusion. The linear 
growth of (a; 2 ) with t, which is characteristic for normal diffusion, results from the Central Limit Theorem applied 
to many independent jumps of a random walker. The anomalous diffusion occurs when the famous theorem fails 
to describe the system because the distributions of summed random variables are too broad or the variables are 
correlated to each other. In physical terms, the subdiffusion is related to infinitely long average time that a random 
walker waits to make a finite jump. Then, its average displacement squared, which is observed in a finite time interval, 
C*~) | is dramatically suppressed. 

The subdiffusion has been recently extensively studied, see e.g. [1-7]. While the phenomenon is theoretically rather 
well understood there are very a few experimental investigations. There is no effective method to experimentally 
measure the parameters a and D a . In the pioneering study [7], where the subdiffusion coefficient was determined 
experimentally for the first time, the interdiffusion of heavy and light water in a porous medium was observed by 
means of nuclear magnetic resonance. The subdiffusion coefficient D a was determined, using the special case a = 2/3 
solution of the fractional derivative diffusion equation. The procedure is neither very accurate nor of general use. 
Our aim here is to develop a method to precisely measure the parameters a and D a . For practical reasons, which 
O ■ are explained below, we choose for the experimental study a membrane system containing two vessels with a thin 
membrane in between which separates the initially homogeneous solute of the substance of interest from the pure 
solvent. A schematic view of the system is presented in Fig. 1. We show that the membrane does not affect, as 
expected, the values of investigated parameters. Instead of the mean square displacement (1), our method refers to 
the temporal evolution of the thickness 8 of the so-called near-membrane layer which is defined as a distance from the 
membrane where the substance concentration drops n times with respect to the membrane surface; k is an arbitrary 
number. In our previous paper [8], we demonstrated that for the normal diffusion S(t) = A\ft. Here we are going to 
study the diffusion of glucose and sucrose in a gel solvent to show that 

5{t) = AV, (2) 

with 7 < 1/2. Our choice of the gel medium to look for a subdiffusion is far not accidental. A gel is built of large 
and heavy molecules which form a polymer network. Thus, the gel water solvent resembles a porous material filled 
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with water. Since a mobility of sugar molecules is highly limited in such a medium the subdiffusion is expected. As 
we show below, this expectation is indeed fulfilled. 

Using an analytic solution of the fractional subdiffusion equation, we theoretically argue, that a = 2-f and that the 
parameter A uniquely depends, as in the case of the normal diffusion, on the number k and the parameters a and D a . 
Therefore, knowing A, n and a, one can deduce D a for arbitrary a. Since the subdiffusion parameter a is measured 
with a high accuracy we can distinguish a very slow normal diffusion from a subdiffusion with a close to unity. Our 
very first theoretical and experimental results on the subdiffusion parameter a have been published in two short notes 
[9]. 

II. EXPERIMENT 

The membrane system under investigation is the cuvette of two vessels separated by the horizontally located 
membrane. Initially, we fill the upper (lower) vessel with the solute of transported substance while in the lower 
(upper) one there is a pure solvent. Then, the substance diffuses from one vessel to another through the membrane. 
Since the concentration gradient is in the vertical direction only, the diffusion is expected to be one-dimensional (along 
the axis x). 

The substance concentration is measured by means of the laser interferometric method. The laser light is split into 
two beams. The first one goes through the membrane system parallclly to the membrane surface while the second 
(reference one) goes directly to the light detecting system. The interferograms, which appear due to the interference 
of the two beams, are controlled by the refraction coefficient of the solute which is turn depends on the substance 
concentration. The analysis of the interferograms allows one to reconstruct the time-dependent concentration profiles 
of the substance transported in the system and to find the time evolution of the near-membrane layers which are of our 
main interest here. We note that the measurement does not disturb the system under study. The experimental set-up 
is described in detail in [10]. Here we only mention that it consists of the cuvette with membrane, the Mach-Zehnder 
interferometer including the He-Ne laser, TV-CCD camera, and the computerized data acquisition system. 

For each measurement, we prepared two gel samples: the pure gel - 1.5% water solution of agerose and the same 
gel dripped by the solute of glucose or sucrose. The concentration of both sugars in the gel was fixed to be either 
0.1 [mol/dm 3 ] or 0.07 [mol/dm 3 ]. The two vessels of the membrane system were then filled with the samples and the 
(slow) processes of the sugar transport across the membrane started. We used an artificial membrane of the thickness 
below 0.1 mm. The membrane was needed for two reasons. It initially separated the homogeneous sugar solute in 
one vessel from the pure gel in another one. It also precisely fixed the geometry of the whole system. 

When the sugar was diffusing across the membrane we were recording the concentration profiles in the vessel which 
initially contained pure gel. The examples of typical interferograms and extracted concentration profiles are presented 
in the earlier paper [10]. The thickness of a near- membrane layer 8 was calculated from the measured concentration 
profiles C(x,t) according to the definition 

C(5,t) = KC(0 + ,t) , (3) 

where k is an arbitrary number smaller than unity. In our analysis we used k = 0.12, 0.08 and 0.05. Computing 6 for 
various t, we found the thickness of near-membrane layer as a function of time. 

In Fig. 2 we present 5 as a function of time for the glucose and sucrose transported in a gel for the initial sugar 
concentration equal 0.1 [mol/dm 3 ]. For the glucose we took three values of k = 0.12, 0.08 and 0.05 while for the 
sucrose k = 0.08. The errors of S(t) shown in the figure were found in the following way. We performed 6 independent 
measurements of the concentration profiles C(x,t) at the same initial sugar concentration. The thickness of the near 
membrane layer 5(t) was found for each concentration profile. We note here that S(t) is insensitive to a sizeable error 
of absolute normalization of C(x, t) as it cancels due to the definition (3). Having 6 values of S for every t, we obtained 
the mean value of S and the standard error. The final errors shown in Figs. 2 were obtained by further multiplying the 
standard errors by the Student-Fisher coefficient taken at a confidence level 95% to include the effect of low statistics. 

As seen in Fig. 2, the experimentally found time dependence of S is well described by the power function (2) with 
the common index 7 = 0.45. In Fig. 3 we show the same data as in Fig. 3 but in a log-log scale to better test the 
power law dependence. The lines representing \ft are also shown for comparison. It is evident that the measured 
index 7 is smaller than 0.5. There are some deviations of our data from the power law dependence at t < 300 s. 
However, as explained in Sec. IIIB, our theoretical formulas, in particular the power law behavior, are derived in the 
long time approximation. 

We fitted the experimental data shown in Figs. 2, 3 by the power law formula At 1 , using a standard \ 2 procedure. 
The values of \ 2 per degree of freedom were smaller than unity for every value of k. The error of the universal index 
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7 = 0.45 was found to be 0.005. The parameter A depends on k; it is also different for glucose than for sucrose. We 
found 



A = 0.091 ± 0.004 for k = 0.05 , (4) 
A = 0.081 ± 0.004 for k = 0.08 , (5) 
,4 = 0.071 ± 0.004 for k = 0.12 , (6) 

for glucose and 

A = 0.064 ±0.003 for k = 0.08 . (7) 

for sucrose. 

The fact that 7 is the same for glucose and sucrose suggests that the index depends only on the solvent i.e. the 
medium where the subdiffusion occurs. As will be shown in Sec. IV, the subdiffusion coefficients for glucose and for 
sucrose differ from each other. Such a situation resembles the case of normal diffusion where 7 = 1 is universal but 
the diffusion coefficient changes from one substance to another. 

At the end of this section we note that the results of measurements of the near-membrane layers for the initial sugar 
concentration equal 0.07 [mol/dm 3 ] fully coincide with those presented above which were performed for the initial 
concentration equal 0.1 [mol/dm 3 ]. This simply reflects the linearity of the subdiffusion equation which is discussed 
in the next section. 

III. THEORY 

In this section we discuss a theoretical description of the experimental result presented in Sec. II. In particular, we 
derive the formula which allows one to obtain the subdiffusion coefficient D a from the experimental data. 

A. General formulas 

The subdiffusion is described by the equation [1,11] 

ac(x,t) _ a 1 -" d 2 c( Xl t) 
dt a d^- a dx* ' [ ' 

with the Riemann-Liouville fractional derivative defined as 

d a .f(t) 1 d n r* f{t') 



- a) dt n J (t - f 



dt a Y{n-a)dt n J (t - f ')!+«-« 
where n is the smallest integer larger than a > 0. It is tempting to write down Eq. (8) in a simpler form as 

d a C(x,t) _ d 2 C(x,t) 

dt« - a dx* ■ iyj 

In general, however, the fractional derivative does not obey the property 

and Eq. (8) is not fully equivalent to Eq. (9) [12]. And because the subdiffusion equation as derived in the framework 
of Continuous Time Random Walk [11] is of the form (8), it has become customary to use Eq. (8) rather than Eq. (9). 
We also note here that Eq. (8) with the temporal fractional derivative of order a < 1 corresponds to a 'long' (infinite 
on average) waiting time of the random walker. This is just the physical situation which is expected in a gel medium 
built of large and heavy molecules forming a polymer network where mobility of the walker is strongly limited. 
We solve Eq. (8) in the region x > for the following initial concentration 

^o) = H° ; x x <° > do) 
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Here, x = is the position of an infinitely thin membrane. In fact, we solve Eq. (8) for the Green's function G(x, t: xq) 
with the initial condition 



G(x, t = 0; x ) = S(x — x ) , (11) 
and then, the concentration profiles are calculated using the integral formula 



C(x, t)= J dx G(x, t; x ) C(x , 0) . (12) 

The Green's function G(x, t; x ) gives the probability density to find a random walker at the position x in time t; the 
walker starts from xo at t = 0. 

To find the concentration profile C(x,t) and then time evolution of near- membrane layer 5, we use the relation 

G(x,t;x )= [ dt' J(0 + ,i';x o )G rcf (x,i-i';0 + ) , (13) 
Jo 

where x > while xo < 0; J(x, t; x ) is the flux associated with G(x, t; xq) which for x = gives the flow across the 
membrane; G ro f(x, t; xq) is the Green's function for the half-space system with x > and the fully reflecting wall at 
x = 0. The integral formula (13), which we call the wall relation, represents one-half of the membrane system as a 
half-space system with a reflecting wall replacing the membrane. The substance flux, however, does not vanish at the 
wall but it equals the actual flux in the membrane system. For the case of normal diffusion the wall relation (13) was 
used in Ref. [14]. In Appendix A we show that it also holds for subdiffusion. 
Using the relation (13) and Eq. (12), the concentration profile can be written as 

C(x,t)= ( dt'W(t')G Icl (x,t-t';0 + ) , (14) 
Jo 

where x > and the function W(t), which equals 

,o 

W(t) = dx J(0 + ,t;x o )C(xo,0) , (15) 

J — oc 

contains all information about the initial and boundary conditions. The interval of integration in Eq. (15) is due to 
the initial condition (10). 

Since the subdiffusion equation is of the second order with respect to the space variable x, it requires two boundary 
conditions at the membrane. The first one simply assumes the continuity of the flux J flowing through the membrane 
i.e. 

J(0-,t) = J(0+,t), (16) 

where J is the subdiffusion flux given by the generalized Fick's law [15]. It gets a simple form after Laplace transfor- 
mation L {/(t)} = f(s) = J °° f(t)e- st . Then, it reads [15] 

J(x, S ) = -^.s 1 -"^^. (17) 

There is no obvious choice of the second boundary condition. Therefore, we assume that the missing condition is 
given by a linear combination of concentrations and flux at the membrane i.e. 

biCiO- ,t) + b 2 C(Q + ,t) + b 3 J(0,t) = 0. (18) 

Since the current is continuous at the membrane there is no point to distinguish J(0 + , t) from J(0~, t). We note that 
two boundary conditions discussed in the literature are of the general form (18). The first condition can be formulated 
as follows: If during a given time interval, N particles reach the membrane, the fraction a of them will be stopped 
while (1 — a) will go through. This condition leads to the boundary condition of the form [8,16-18] 

C(0+,t) = I^<7(0-,t). (19) 
1 + er 
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The second condition assumes that the flux flowing through the membrane is proportional to the difference of con- 
centration at membrane surfaces [19,20] 



J(0,t) = X(C{Q-,t) - C(0+,i)) 



(20) 



where the parameter A controls the membrane permeability. We call the boundary conditions (16,19) and (16,20) as 
'A' and 'B', respectively, and in Sec. IIIC, we use the indices 'A' and 'B' to denote the solutions of the fractional 
subdiffusion equation obeying the corresponding boundary conditions. 

The initial condition (10) combined with the general boundary condition (18) give the Laplace transform of the 
function W(t), which is defined by Eq. (15), as 



W(s) = Co 



h-b 2 - hy/D^s 1 -"/ 2 W 2 



I =Co W^_l_g dfesfe (l-a/2) 



where 



'b 1 -b 2 W 2 



'D n . 



(21) 



k=0 



bi - b 2 

The derivation of Eq. (21) is discussed in Appendix B. Inverting the Laplace transform, we get 

b iy /D^ 1 A d k 1 



W(t) = C Q 



fl-a/2 



h - b 2 i 1 -"/ 2 ^ r(a/2 - fc(l - a/2)) i*(i-«/2) ' 



(22) 



The Green's function G re f, which enters Eq. (14), can be easily obtained by means of the method of images [1] as 

G lei (x, t; x ) = G (x, t; x ) + G (-x, t; x ) , (23) 
with Go being the Green's function for the homogeneous system given by 

1JD 



G (x,t;xo) 



a\x — xo 



Tjl 
H l 1 



o f f \x -xp\ 



1 1 

1 1 



where H denotes the Fox function. In our numerical calculations we use the form of H 

/ i \ 1 + lJ 



1 o 



rrl U 

H l 1 



a" 



1 1 a 1 a " 

l±v 1 = P 

a 



t 



f^ T(-kp - v)k\ V tP) 



(24) 



(25) 



which is derived in [13] by means of the Mellin transform technique [22]. In terms of the Laplace transform, the 
function (24) simplifies to the form [1] 



Go(x, s;x ) 
which was found using the formula [13] 

L-\s v e- a,fi ) 



\x— Xo|l 



(3a- 



1 rrl Of 0* 



1 1 

l+v 1 





(26) 



(27) 



where a, (3 > 0. 

Having the explicit form of the functions W and G rc f , we can write down the concentration profile using Eq. (14). 
The Laplace transform of G equals 



d(x,s) -Co- 



61 



1 



bi -b 2 - b 3 y/D^ s 1 -*/ 2 s 
which after inverting the transformation can be written down as 



C(x, t)= [ dt' W(t - t') — Hi ° 
Jo ax 



1 1 



(28) 



(29) 



Eq. (29) is the starting point of our further analysis which provides the function 5(t). 
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B. Long time approximation 



We first consider the long time approximation of the formula (29) which, according to the Tauberian theorem [21], 
corresponds to the small s limit of the Laplace transform. The physical meaning of this approximation is discussed 
below. Taking into account only the leading contribution in the small s limit, Eq. (28) gets the form 



C{x,s)=C - , 

bi - 6 2 s 

which after inverting Laplace transformation provides 



C(x, t) = Co 



2&i 



(h - b 2 )a 



ttI o 

H i l 




(30) 



(31) 



The solution (31) can be also obtained directly from Eq. (29), taking into account only the k = term in the expansion 
of W{t) given by Eq. (22), and using the formula 



dt' 



ttI 
H l 1 



xT(a/2)J {t-t'Y 



-a/2 



i i 



1 



rrl U 

= H l 1 




(32) 



which is derived by means of the Laplace transformation. 

The series (22) can be approximated by the first term if d <C t 1 ^"^ 2 which gives 



< t . 




(33) 



When the boundary condition is of the form (19), the condition (33) is trivially satisfied for any t as 63 = in this 
case. Even more, as discussed in the next section, the solution (31) is exact for the boundary condition (19). When 
the boundary condition is of the form (20), we have A = 61/63 = —62/63, and the condition (33) can be written as 



2A 



— 



< t . 



(34) 



We have first estimated the l.h.s. of Eq. (34), assuming that we deal with the normal diffusion (a = 1). For the 
membranes used in our experiments the parameter A is of order 10~ 2 [mm/s] while the coefficient of normal diffusion 
D is roughly 10~ 5 [mm 2 /s]. Thus, the l.h.s. of Eq. (34) is estimated as 2 s. Since 10 s is the time step of our 
measurements which extend to 2500 s, the condition (34) is fulfilled. We have also checked the condition (34) a 
posteriori, using the values of a and D a found in Sec. IV. The l.h.s. of Eq. (34) is again about 2 s. 

Let us now discuss the temporal evolution of near-membrane layers in the long time approximation. Substituting 
the solution (31) into Eq. (3), which defines the near-membrane layer, we get the equation 



1 



H l 1 




= kH{ ? 1 



1 1 

* 



(35) 



which due to the identity 



simplifies to 



1 1 

# 



Q 

2 



Tjl 
H l 1 



One observes that Eq. (36) is solved by 




6(t) =A(a,D a , K ) t a ' 2 



(36) 



(37) 
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where the coefficient A equals 



A(a, D ai k) 



(Hi",)- 1 



an 
~2~ 



1 1 

# 



(38) 



We note that the near-membrane layer given by Eq. (37) does not depend on parameters b\ and 62 which control 
membrane permeability. 5(t) remains the same in the absence of membrane. We also observe that the coefficient A 
can be recalculated into the diffusion constant D n as 



A 2 



1 f 

I 



(39) 



Thus, knowing experimental values of A, a and k, one can deduce D a 



C. Specific boundary conditions 

In the previous section, we have determined the temporal evolution of near-membrane layer for arbitrary linear 
boundary condition (18) in the long time approximation. Here we discuss the problem for two specific sets of the 
boundary conditions (16,19) and (16,20) which are called 'A' and 'B', respectively. The solutions of the subdiffusion 
equation (8) satisfying the boundary conditions A and B can be obtained from the general solution (28), substituting 
b\ = a + 1, 62 = c — 1, 63 = or &i = A, b 2 = —A, 63 = —1, respectively. However, the solutions of Eq. (8) with 
the boundary conditions A and B can be also found in a much simpler as demonstrated in Appendix C. 

As in the case of normal diffusion, the boundary conditions (16,19) allow for analytic solution of Eq. (8). Using the 
technique of Laplace transform (see Appendix C for details), the concentration profile for x > is found as 



C A {x,t) = {l-a 




1 1 




2 I • (40) 



The solution (40) plugged into the definition of near-membrane layer (3) gives Eq. (36) which is solved by S\(t) of 
the form (37). Thus, as already mentioned, the formulas derived in the long time approximation are exact for the 
boundary condition A. 

When Eqs. (16, 20) are used as the boundary conditions, the solution of subdiffusion equation (8) can be written 
in the form of the infinite series of Fox functions. As explained in Appendix C, for x > one finds 



„ 00 



x 
2A 



rrl 

n=0 



D a t c 




(41) 



The solutions (40) and (41), which for normal diffusion have been discussed in [16-19], qualitatively differ from 
each other. In particular, according to Eq. (40) the flux flowing through the membrane is constant in time while the 
solution (41) leads to the flow which decreases in time as the concentrations at both sides of the membrane evolve to 
the same value Co/2. In the case of Eq. (40), the ratio of concentrations at both sides is constant in time as dictated 
by the boundary condition A. However, the qualitative differences of the solutions (40) and (41) are evident for times 
which are so long that substance concentration becomes approximately homogeneous at each side of the membrane 
in its neighborhood. These times are usually much longer than those satisfying the condition (33) when the vessels 
of a membrane system are of a few centimeters length as those used in our measurements. 

Since the formula (41) is analytically intractable we have found the time evolution of near-membrane layer numeri- 
cally. The results are shown in Figs. 4-8. The physical meaning of the near-membrane layer suggests that it does not 
depend on the membrane permeability. The analytical calculations performed with the boundary condition A fully 
confirm this intuition. As seen, 5a remains the same even in the absence of the membrane (<r = 0). We are now going 
to show numerically that the same holds for the boundary condition B. Specifically, we argue that <5a(£) = #b(£) for 
the same values of k, a and D a , in spite of the qualitative differences between the solutions (40) and (41). 

Fig. 4 demonstrates that <$b(£) is independent of the membrane permeability parameter. Our numerical calculations, 
which are presented in Figs. 5-8, show that the thickness of the near-membrane layer grows in time as V with 7 = a/2. 
Thus, we write 

S B (t)=A B t~<. (42) 
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We are now going to show that Aa = Ab- In Figs. 5-7, we compare $B(t) obtained numerically with #a(£) given by 
Eq. (37) for the same values of k, a and D a . In each of these figures one parameter changes whereas the remaining 
two are fixed. In the Fig. 5 we examine the dependence on the parameter defining thickness of near membrane layer 
k; in Fig. 6 there are several values of the subdiffusion coefficient D a while in Fig. 7 the functions are plotted for 
several values of a. In all these cases we see a perfect coincidence of 5 a (t) with 5b (t). Fig. 8 demonstrates that 
Ssit) fulfills the scaling exactly obeyed by 5 a (t). Namely, we show that 5'(t) = Sb^/Aa with A a given by Eq. (38) 
depends solely on time, exactly as 5a(£)A4a- 

We conclude this section by saying that no difference between 5a (t) and <5b(0 has been observed. Therefore, 
we expect that Aa = Ab = A. Consequently, the subdiffusion coefficients calculated from Eq. (39) have to agree 
with the ones obtained numerically form the solution (41). Although we are unable to prove beyond the long time 
approximation that the time evolution of near-membrane layers does not depend on the boundary condition, the 
results presented in this section strongly substantiate such a conjecture, and thus they justify the use of the formula 
(39) to evaluate the subdiffusivc coefficient D a from experimental data. 



IV. EXPERIMENTAL VALUES OF a AND D a 

In Sec. II we have fitted the experimentally obtained 5(t) by the power function At 1 . Thus, we have found the 
index a — 2j — 0.90 ± 0.01 and the values of the coefficient A given in Eqs. (4,5,6) and (7). Now, we recalculate A 
into D a by means of the relation (39). Using the numerical values of inverse Fox functions 



0.054 



0.036 



1 1 

ois 



0.0225 



1 
i 

0.45 

1 

1 

0.45 



= 6.032 



8.014 



10.510 



which are calculated with the help of the expansion formula (25), we computed the subdiffusion coefficient for each 
value of k. As expected, D a shows no dependence on k within the errors which were obtained propagating the errors 
of A. We took the mean value of D a as a final result but the final error is the maximal one. The point is that the 
values of D a found for different values of k cannot be treated as independent measurements. Thus, we obtained 



for glucose and 



Ai.go = (9.8 ± 1.0) x 10~ 4 [mm7s - 90 ] 



D .9o = (6.3 ± 0.9) x 10~ 4 [mm 2 /s° 



for sucrose. 

To be sure that Eq. (38), which is used to evaluate D a , properly describes the experimentally found 5(t), we 
checked the scaling property of 5(t) suggested by the formula (37). In Fig. 9 we plotted the rescaled thickness of 
near-membrane layer S'(t) = 5(t)/A with A given by Eq. (38) for all values of k, for glucose and for sucrose. The 
experimental points are represented as in Fig. 2. According to Eqs. (37,38), 5'(t) is simply the function t 7 , and, as 
seen in Fig. 9, all our experimental data are indeed very well described by t - 45 . 



V. FINAL REMARKS 



Our method to determine the parameters of subdiffusion relies on the near-membrane layers. One may ask why a 
and D a arc not extracted directly form the concentration profiles which are measured. There are three reasons to 
choose the near-membrane layers: experimental, theoretical and practical. 

1. Measurement of the near-membrane-layer thickness does not suffer from the sizable (~ 10-15%) systematic 
uncertainty of absolute normalization of the concentration profiles, as only the relative concentration matters 
for (5, see the definition (3). 
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2. Computed concentration profiles depend on the adopted boundary condition at a membrane while the condition 
is not well established even for the normal diffusion. The near-membrane layer is argued to be free of this 
dependence. 

3. When the concentration profile is fitted by a solution of the subdiffusion equation, there are three free parameters: 
a, D a and the parameter characterizing the membrane permeability. Because these fit parameters are correlated 
with each other it is very difficult in practice to get their unique values. One should remember here that the 
solution of the subdiffusion equation is of rather complicated structure, see Eqs. (40, 41), which makes the 
fitting procedure very tedious. When the temporal evolution of 8 is discussed the membrane parameter drops 
out entirely, a is controlled by the time dependence of S(t) while D a is provided by the coefficient A. 

We have argued using Eq. (8) that 8 scales as V . A similar scaling can be obtained in diffusion equation with 
fractional derivatives in space and time if the orders of the fractional operators are properly chosen. This is related to 
the fact that in fractional diffusion equation of order /i in space and order v in time, (x 2 ) scales as t 2l/ > ^ . One could 
also come up with a nonlinear diffusion model without fractional derivatives that exhibits the observed non-Gaussian 
scaling for 8 [23]. This non- uniqueness problem cannot be solved in full generality but some possibilities can be 
excluded. The fractional differential equation of the order \x < 2 of spatial derivative corresponds to 'long' jumps (the 
variance is infinite) of the random walker. Since we study experimentally the diffusion in a gel built of large and heavy 
molecules forming a polymer network, long jumps are physically not plausible. We rather expect that the average 
waiting time of the random walker in a gel medium goes to infinity. This corresponds to the fractional differential 
equation with the temporal derivative of order a < 1 as in Eq. (8). One cannot exclude the equation with fractional 
spatial and temporal derivatives. Similarly, the nonlinear diffusion equation cannot be excluded. However, Eq. (8), 
which is derived in a physically well motivated Continuous Time Random Walk approach [1,11], offers the simplest 
possibility. We also note that the scaling shown in Fig. 9 strongly supports the subdiffusive model based on Eq. (8). 
It seems difficult to reproduce this scaling in a quite different theory. 

In our calculations presented in Sec III, the membrane is assumed to be infinitely thin. Relaxing this assumption 
considerably complicates theoretical analysis of the problem. This is not only the finite membrane thickness that 
should be taken into account but the membrane internal structure should be modeled as well. In particular, one 
should answer the question whether transported substance is accumulated inside the membrane. If so, the membrane 
permeability might be time dependent. To avoid all these complications, the membrane is infinitely thin in our analysis 
which is a reasonable assumption as long as the membrane is sufficiently thin. (The thickness of the membrane used 
in our measurements was below 0.1 mm.) However, we expect that our method to determine the parameters a and 
D a also works for finite width membranes. First of all, the temporal evolution of near-membrane layer was shown 
to be fully independent of the boundary condition at a membrane in the long time approximation. Our numerical 
calculations also suggest that 8{t) does not depend on the membrane. Finally, we observe that for the near-membrane- 
layer (3) only the relative substance concentration at x — 8 and at the membrane surface (x — 0) matters. Therefore, 
it is expected that the membrane properties do not influence the time evolution of near-membrane layer. 

Our method to determine the subdiffusion parameters uses the membrane system. While the membrane plays here 
only an auxiliary role, it should be stressed that the transport in membrane systems is of interest in several fields of 
technology [24], where the membranes are used as filters, and biophysics [25], where the membrane transport plays 
a crucial role in the cell physiology. The diffusion in a membrane system is also interesting by itself as a nontrivial 
stochastic problem, see e.g. [16-18,20]. Thus, our study of the subdiffusion in a membrane system, which to our best 
knowledge has not been investigated by other authors, opens up a new field of interdisciplinary research. It is also 
worth mentioning that our intcrfcrometric set-up can be used to experimentally study the (anomalous) diffusion not 
only in the membrane systems. In particular, we plan to perform measurements in a system with no membrane where 
the sugar is transported directly from the water to gel solvent. However, there are problems to keep fixed geometry 
of the whole system in the course of long lasting diffusion process. 

At the end let us summarize our considerations. We have developed a method to extract the subdiffusion parameters 
a and D a from experimental data. The method uses the membrane system, where the transported substance diffuses 
from one vessel to another, and it relies on the fully analytic solution of the fractional subdiffusion equation. We have 
applied the method to our experimental data on the glucose and sucrose subdiffusion in a gel solvent and we have 
precisely determined the parameter a and the subdiffusion coefficient D a . 
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APPENDIX A: WALL RELATION 



We derive here the integral relation (13) which after the Laplace transformation reads 

G+-(x,s;x ) = J + _(0 + ,s;x o ) G rof (x, s; 0+) . (Al) 

We have introduced here the indices + and — which correspond to the signs of x and x . Since x > and xo < 0, 
the Green's function is labeled with H — . 

Because the current J(x, s;x ) is expressed, in accordance with Eq. (17), as 

J(x,s;x ) = -D aS ^ , (A2) 

and the general solution of the subdiffusion equation is, as discussed in the Appendix C, given by Eq. (C2), the current 
equals 

J+_(0+ s;a;o) = yfD^s 1 - a/2 B{s) + X - e x °V^ = /d^ 1 -"/ 2 G+_(0+, s; x ) . (A3) 
As the Green's function of the system with the reflecting wall at x = is given by Eqs. (23,24,26), we have 

s; 0+) = vz^ 1 -"/ 2 e ~°^ ' (A4) 

and due to Eq. (A3) the expression in the r.h.s. of Eq. (Al) equals 

J + _(0 + ,s;x o )G rcf (a;,s;0 + ) ^e^V^ G + _ (0+, s; x ) . (A5) 
Eq. (C2) allows one to write down the Green's function from the l.h.s. of Eq. (Al) as 

G+-(x, S ;x ) = e^VSF G + _ ( +, 3 - Xq ) . (A6) 
Comparing Eq. (A5) to Eq. (A6), one finally finds the wall relation (Al). 



APPENDIX B: FUNCTION W 



We derive here Eq. (21). Taking the Laplace transform of the boundary condition 18, replacing the concentration 
profiles by the Green's functions and putting xq < 0, we get 



b 1 G—{0-,s;x ) + b2G + -{0-,s;x ) + b 3 J(0,s;x ) = 



Using Eq. (C2), one finds 



G (x,s;xq) = A(s) e x V D ° 



G + _ (x, s-xq) = B(s) e 



1 



1 



+ ■ 



-\x— Xo|l 



-\x— Xo|a 



2y/D a s 1 -^ 

The Green's functions (B2,B3) provide, via Eq. (17), the currents 

1 
2 
1 



J — (0 ,s;xq) 



D a A{s) 



J + _(0 + ,s;x o ) = -e X0 V n a + ^ Da B(s) . 



(Bl) 

(B2) 
(B3) 

(B4) 
(B5) 



Due to the current continuity equation (16), A(s) = —B(s) and the current, which enters the boundary condition (Bl), 

equals J H (0 + , s;x ). Substituting the expressions (B2, B3, B5) into Eq. (Bl) and using the equality A(s) = —B(s), 

one finds 



B(s) = 



b 1 +b 2 + b z ^D2s l - a l 2 



(B6) 



2VA^s 1 -«/ 2 (6 1 - b 2 - feaV^s 1 -"/ 2 ) 

The function B(s) determines the current (B5) which, substituted into the definition (15) together with the initial 
condition (10), provides the function W(s), as given by Eq. (21). 
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APPENDIX C: SOLVING SUBDIFFUSION EQUATION 



We briefly present here a procedure of solving the subdiffusion equation (8) with the boundary conditions (16, 19) 
and, respectively, (16, 20). Taking the Laplace transform of Eq. (8), we obtain 



C(x, s) - s- a D a d2C ^ S) ~ C(x, 0) = , 



(CI) 



where, as previously, we use the hats to denote the Laplace transformed functions. The solution of Eq. (CI) with the 
initial condition (11) reads 



G(x, s; x ) = A(s) e x V D ° + B(s) e X V D = 



1 



-\x-xo\W -fe- 



2 V / I^s 1 ^ 

where the functions A(s) and B(s) are determined by the boundary conditions. 



(C2) 



1. Boundary conditions A 



The Laplace transform of boundary condition (19) is 



C(0 + ,s) = ]—^C(0-,s) 
1 + a 



(C3) 



Proceeding analogously to the case of normal diffusion [17,18], we find the Green's functions obeying the boundary 
condition (C3) as 



G — (a;,s;a;o) = G Q (x, s;x a ) + aG (-x, s; x ) , 
G+-(x, s;x ) = (1 - a) G Q (x, s;x ), 



(C4) 



where Go is the Laplace transform of Green's function for the homogeneous system without a membrane (24); the 
indices + and — of the Green's functions refer to the sign of x and of xq , respectively. 

To compute the concentration profiles we use the Laplace transformed of the integral relation (12) which takes the 
form 



C(x,s)= J G(x,s;x ) C(xo,0)dxo ■ 
Now, we substitute the Green's function (C4) into Eq. (C5) and use the formula 



"T+77 #n — 



(3a 



1 1 

l+v 1 



(C5) 



(C6) 



which is derived in [13] using the Mellin transform technique [22]. Here, a, (3 > while the parameter v is not limited. 
After simple calculations, we finally find the solution (40). 



2. Boundary conditions B 

The subdiffusion equation (8) is solved by the Green's functions of the form 
G — (x, s; xq) = Ai(s) e x V^ + B 1 (s) e~ x ^ + 



1 



G+-(x,s;x Q ) = A 2 (s)e x V D « + B 2 {s) e 



+ 



1 



-\x— Xo|i 



-|x— Xo|l 



(C7) 
(C8) 



The functions A\(s), A 2 (s), B\(s), and B 2 (s) are now determined by the boundary conditions (16,20) which after the 
Laplace transformation take the form 
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J(0 ,s;x ) = J(0 + ,s;x o ) , 
7(0, s; so) = A (6—(0-,s;x ) - G+_(0 + , a; sc )) 
with the Laplace transform of subdiffusive flux given by the formula [15] 



(C9) 
(CIO) 



j(x, s; xq) = —D a s 



a dG(x,s;x ) 
dx 



For the infinite system, one demands vanishing of the Green's functions for x — ► ±00, which gives B\ = A2 = 0. 
Substituting the Green's functions (C7, C8) into the boundary conditions (C9, CIO), we obtain 



G — (x, s;x ) = 
G+-(x,s;x ) = 



4A + 2VA^s 1 -* 
4A 



1 



1 



-(x-x )a/#- 



-\x-x \W -5^ 



4X + 2VD^s 1 -- 2VD^s 1 ~2 



(Cll) 



Expanding the function (Cll) into the power series with respect to s and using the initial condition (10), the integral 
relation (C5) provides, with the help of the formula (C6), the solution (41). 
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FIG. 1. Schematic view of the membrane system under study. 
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FIG. 2. The experimentally measured thickness of the near- membrane layer 5 as a function of time t for glucose with k = 0.05 
(□), k — 0.08 (o), k = 0.12 (A), and for sucrose with k = 0.08 (0). The lines represent the power function At 0A5 with the 
coefficients A given by Eqs. (4-7). 
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FIG. 3. The same experimental data as in Fig. 2 but in a log-log scale. The solid lines represent the power function A i ' 45 
with the coefficients A given by Eqs. (4-7). The dotted lines correspond to the function A\ft. 
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FIG. 4. The time evolution of the near-membrane layer 8b obtained from Eq. (41) for different values of the membrane 
permeability parameter A: 1 [mm/s] (solid line), 5 [mm/s] (□), 100 [mm/s] (x), and 500 [mm/s] (0). The remaining parameters 
equal k = 0.12, a = 0.90, and Do.go = 3 x 10~ 3 mm 2 /s ' 90 . 
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FIG. 5. Comparison of 5a (t) (symbols) and &(£) (lines) 
(V)- The remaining parameters equal a = 0.90 and Do.90 ; 



for several values of k: 0.05 (□), 0.08 (o), 0.12 (A), 0.20 (0), 0.50 
= 3x 10- 3 mm 2 /s ' 90 . 
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FIG. 7. Comparison of 6 A {t) (symbols) and 8 B (t) (lines) for several values of a: 0.9 (□), 0.6 (o), 0.4 (A), 0.2 (0). The 
remaining parameters equal k = 0.12 and D a = 3 x 10~ 3 mm 2 /s Q1 . 
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FIG. 8. The thickness of the near-membrane layer 5b and the rescaled thickness 5b /Aa as functions of time for several values 
of k: 0.05 (□), 0.08 (o), 0.12 (A), 0.20 (0). For the rescaled thickness the symbols are on top of each other. The parameters 
a and D a equal, respectively, 0.90 and 3 x 10 -3 mm 2 /s a . 




FIG. 9. The experimentally measured 8 divided by the coefficient A given by Eq. (38). The symbols are assigned as in Fig. 2. 
The parameters equal a = 0.90 and D .go = 1.0 x 10~ 3 [mm 2 /s ' 90 ] for glucose and D . 90 = 6.3 x 10~ 3 [mm 2 /s 0,90 ] for sucrose. 
For clarity of the plot the error bars are not shown. The line represents the function i 0,45 . 
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